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ABSTRACT 

Several parametrizations for overshooting in ID stellar evolution calculations coexist in the literature. These parametrizations are 
used somewhat arbitrarily in stellar evolution codes, based on what works best for a given problem or even for the historical reasons 
related to the development of each code. We point out that these different parametrizations correspond to different physical regimes of 
overshooting, depending on whether the effects of radiation are dominant, marginal, or negligible. Our analysis is based on previously 
published theoretical results, as well as on multidimensional hydrodynamical simulations of stellar convection where the interaction 
between the convective region and a stably stratified region is observed. Although the underlying hydrodynamical processes are the 
same, the outcome of the overshooting process is profoundly affected by radiative effects. Using a simple picture of the scales involved 
in the overshooting process, we show how three regimes are obtained, depending on the importance of radiative effects. These three 
regimes correspond to the different behaviors observed in hydrodynamical simulations so far and to the three types of parametrizations 
used in ID codes. We suggest that the existing parametrizations for overshooting should coexist in ID stellar evolution codes and 
should be applied consistently at convective boundaries depending on the local physical conditions. 
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1. Introduction 

For more than 50 years, convection has been described in ID 
stellar evolution codes by the mixing-length theory (MET, cf. 
Bohm-Vitense 1958). It is well known that one of the short¬ 
comings of MET is its inability, by construction, to describe the 
boundary layer between the convective region and a neighboring 
stably stratified region. In the MET picture, the flow simply stops 
at the boundary. This contradicts basic physics, because inertia 
allows the flow to penetrate the stably stratified region, inducing 
mixing beyond the boundary of the convective region (see, e.g., 
Arnett et al. 2015). 

Such an extra mixing at convective boundaries is routinely 
added in stellar evolution calculations, since it is required to 
reproduce well-established observational features across the 
Hertzsprung-Russel diagram (e.g., Maeder 1975; Matraka et al. 
1982; Schroder et al. 1997; Herwig 2005; Pace et al. 2012; Mon- 
talban et al. 2013). Currently, this extra mixing is parametrized 
in very crude ways in stellar evolution codes, and it relies on free 
parameters that need to be calibrated. This situation is very sim¬ 
ilar to the one of convection itself, with MET relying on the mix¬ 
ing length parameter. In fact, overshooting is one aspect of con¬ 
vection, and a better description of overshooting should eventu¬ 
ally rely on a better description of turbulent convection. A lot of 
effort has been devoted to the development of a non-local theory 
of convection, which, in principle, should allow for a consistent 
description of overshooting (Gough 1977; Stellingwerf 1982; 
Xiong 1986; Kuhfuss 1986; Canuto & Mazzitelli 1991; Canuto 
1992; Gehmeyr & Winkler 1992; Wuchterl & Feuchtinger 1998; 
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Deng et al. 2006). Usually, this is done by looking for a working 
closure in the framework of Reynolds-Averaged Navier-Stokes 
(RANS) equations (see e.g. Canuto (1997); Xiong et al. (1997); 
Canuto & Dubovikov (1998); Deng et al. (2006); Canuto (201 1) 
and references therein). 

Although RANS is the appropriate mathematical framework 
for improving the description of hydrodynamical processes in 
stellar evolution, one should abandon the idea of finding the “ul¬ 
timate” parametrization that would describe overshooting cor¬ 
rectly in any physical situations. Although the word “overshoot¬ 
ing” is widely used to describe extra mixing at convective bound¬ 
aries (and we adopt this terminology as well), it is misleading 
because the mixing cannot be seen simply as resulting from a 
unique, basic, physical process. Mixing at convective boundaries 
likely results from the interplay of several physical processes 
(shear instabilities, g-modes, etc.), the usual picture of plume 
penetration being only one aspect of the problem. A better de¬ 
scription of convective boundary mixing requires first of all a 
better understanding of the nature and efficiency of the mixing 
processes that are taking place in the boundary layer and how 
they are affected by the local physical conditions. This work 
takes a first step in this direction by looking specifically at the 
effect of thermal diffusion. 

In Sect. 2, we describe the three types of parametrization for 
overshooting that are used in stellar evolution calculations. Cur¬ 
rently, these parametrizations are used somewhat arbitrarily. In 
Sect. 3, we develop a physical picture for the effects of radiation 
on the overshooting process that leads us to define three distinc¬ 
tive regimes of overshooting. Hydrodynamical simulations sup¬ 
port our picture and are discussed. Furthermore, we show that 
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the different parametrizations that are used currently correspond 
to these different thermal regimes of overshooting. We conclude 
in Sect. 4 that the description of overshooting in stellar evolution 
codes could be already improved by using existing parametriza¬ 
tions in a more physically consistent way. In case such a better 
practice leads to theoretical predictions contradicting observa¬ 
tions for some specific cases, this would point to shortcomings 
in the current descriptions and provide valuable information for 
improving them. 


2. Existing parametrizations of overshooting 

2.1. “Chemical mixing" 

In this approach, it is assumed that overshooting does not affect 
the thermal structure. In terms of the logarithmic temperature 
gradients, one has V = V^d in the overshooting layer. Regard¬ 
ing the mixing of chemical elements, two approaches are pos¬ 
sible. In the simplest approach, the composition is assumed to 
be mixed instantaneously over a distance lov beyond the limit of 
the convective region, with Iqv typically parametrized in terms of 
the pressure scale height. This is the approach taken in BaSTI 
(Pietrinferni et al. 2004), TGEC (Theado et al. 2012), and DSEP 
(Dotter et al. 2008). Another possibility is to treat this mix¬ 
ing as a diffusive process with a prescribed diffusion coefficient 
Dov Originally, this approach was suggested by Ereytag et al. 
(1996) based on radiative hydrodynamic al simulations of near¬ 
surface convection. Eollowing Ereytag et al. (1996), Dqv is usu¬ 
ally parametrized as 


Doy 


- Do exp 



r-rp 

.fHp 


( 1 ) 


with rp the location of the convective boundary, which is typi¬ 
cally given by the Schwarzschild or Ledoux criteria. Dp, which 
is a reference diffusivity (typically the mixing-length diffusivity 
T^mlt evaluated close to the boundary'), and / a free parame¬ 
ter. This “diffusive mixing” approach is the one taken in MESA 
(Paxton et al. 2011) and GARSTEC (Weiss & Schlattl 2008). 


2.2. “Penetrative convection” 

Zahn (1991) showed that when the convective boundary is lo¬ 
cated deep inside the star, the overshooting material is able 
to change the entropy stratification and induce a nearly adia¬ 
batic (yet subadiabatic) region (see also Schmitt et al. 1984). 
Zahn (1991) coined this “penetrative convection”. In stellar evo¬ 
lution codes, this is implemented by artificially increasing the 
size of the adiabatic region on a distance Igy, which is typically 
parametrized in terms of the pressure scale height^. The major 
difference from chemical mixing is that the thermal structure of 
the star is modified: V = Vad in the overshooting region. In addi¬ 
tion, the chemical composition is assumed to be mixed instanta¬ 
neously in this region. This treatment of overshooting is the one 
implemented in GENEC (Eggenberger et al. 2008). 


' By construction Dmlt is zero at r = rg, so that in practice it is evalu¬ 
ated “close” to the boundary. The definition of Dg likely varies between 
different stellar evolution codes. 

^ Some parametrizations use a maximum limit to prevent problems due 
to the increasing value of Hp toward the center. 


2.3. “Turbulent entrainment” 

Turbulent entrainment is a process that is well-known in geo¬ 
physics, because it is observed both in the atmosphere and in 
the oceans. Entrainment is the physical process by which tur¬ 
bulent eddies entrain mass at the convective boundary, inducing 
a steady growth in the size of the convective region as long as 
energy is supplied to the system. Turbulent entrainment is ob¬ 
served in the simulations of the oxygen-burning shell presented 
in Meakin & Arnett (2007), showing that the relevance of this 
process for stellar evolution has been overlooked. To our knowl¬ 
edge, only Staritsin (2013, 2014) has published stellar evolution 
models that include turbulent entrainment at convective bound¬ 
aries. His approach follows “bulk” entrainment models in which 
the boundary layer is collapsed to a discontinuity separating the 
well-mixed region from the stably stratified region (see, e.g., Eer- 
nando 1991). The velocity Ve at which this boundary moves as a 
result of entrainment can be parametrized as 


— = ARi7", 


( 2 ) 


where Vt is the typical turbulent velocity at the boundary, A and 
n are parameters that characterize the entrainment (see below), 
and Ri^ is the so-called bulk Richardson number, which charac¬ 
terizes the “stiffness” of the boundary. It is defined as 


Rifi - 


lAb 


(3) 


where I is the typical size of the turbulent eddies doing the en¬ 
trainment, and 


= J N^dr (4) 

is the buoyancy jump across the interface (N is the Brunt-Vaisala 
frequency). The spatial integration is performed on a region that 
contains the convective boundary. Once Ve is determined, the 
size by which the convective region is extended because of en¬ 
trainment during one time step is given hy d - VgAt. 

The entrainment law (2) is well-established in fluid dynamics 
(see, e.g., Eemando 1991). The values of A and n were measured 
in different experimental setups and geophysical environments. 
Meakin & Arnett (2007) found that their simulations correspond 
to n ~ 1.05 and A ~ 0.027, which are the values used in Staritsin 
(2013, 2014). 


3. Impact of radiation on overshooting 

3.1. Measure of radiative effects in stellar hydrodynamics 

A characteristic of stellar hydrodynamics is that, at the typical 
temperature and density of stellar plasma, photons are very ef¬ 
ficient in transporting heat. A measure of the effect of radiation 
on the flow is given by the Peclet number Pe, defined as 


time scale for radiative transport of heat 

Pe = ^- — - - -^-——. (5) 

time scale for advective transport of heat 

When Pe » 1, radiation has a negligible impact and the flow 
can be considered as evolving adiabatically. When Pe 1, radi¬ 
ation dominates. In the optically thick interior, heat transport by 
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Fig. 1. Sketch of the different scales of eddies present in the flow, rj 
and A are the Kolmogorov and integral scales, respectively. The scales 
contributing to overshooting are those between and Z,„ax- Depending 
on the location of the transition scale Ij where Pe ~ 1, one can define 
three different regimes of overshooting (see text). 


radiation can be described as a diffusion process, and the corre¬ 
sponding thermal diffusivity (units; cm^/s) is given by 


l6o-T^ 

2>Kp^Cp ’ 


(6) 


where T is the temperature, p the density, k the Rossland opacity, 
Cp the heat capacity at constant pressure (see, e.g., Kippenhahn 
et al. 2012). In this case, the Peclet number can be defined as 


a’ 

where u (resp. 1) is a typical velocity (resp. length) scale of the 
flow. The Prandtl number, defined as the ratio of viscosity to ther¬ 
mal diffusion Vr-vly, takes very low values in stellar plasma, 
typically Pr=10 ® - 10 The situation is very different in geo¬ 
physical flows, where Prandtl numbers are of order unity. In fact, 
many insights from geophysical studies can be generalized to the 
stellar case by introducing the effects of thermal diffusion. 

3.2. Scales involved in the overshooting process 

Stellar convection is characterized by very high values of the 
Reynolds number, implying that the flow is highly turbulent; mo¬ 
tion extends over a wide range of scales (see, e.g., Arnett et al. 
2014). An important characteristic of turbulent convection is that 
transport properties in the bulk of the convective region are dom¬ 
inated by large-scale, coherent structures usually described as 
“convective plumes” (e.g., Stein & Nordlund 1989; Cattaneo 
et af 1991; Porter (&; Woodward 2000; Murphy & Meakin 2011; 
Viallet et al. 2013). 

As convective plumes approach the boundary of the con¬ 
vective region, they are deflected sideways. At the same time, 
inertia allows a plume to penetrate the stably stratified region, 
where it undergoes “buoyancy braking”, which is the process 
in which kinetic energy is converted into potential energy due 
to the work done against gravity in the stably stratified region. 
Although buoyancy braking reduces the kinetic energy of the 
overshooting material, the flow comes to a rest mostly due to the 
turbulent dissipation that arises as it breaks apart and forms a tur¬ 
bulent cascade through pressure-strain effects resulting from the 
interaction with the stable region. We coin this process “buoy¬ 
ancy breaking”. 

If we adopt the picture of a turbulent flow extending from the 
integral scale A down to the Kolmogorov dissipation scale 77 , it 
is sound to assume that overshooting is effectively achieved by a 
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Fig. 2. Radial profiles of the length scales Lo,L^, and the Peclet number 
in the boundary layer of the red giant model of Viallet et al. (2013). The 
thick dashed lines denote the extent of the boundary layer, defined as 
the region of buoyancy braking. The jc-axis is normalized by the radius 
of the star R = 4.1 x 10'^ cm. 


certain subrange of scales comprised between /min and /max (see 
Fig. 1). Here, /max is the typical length scale of plumes, while 
/min is the length scale of the smallest eddies that are able to con¬ 
tribute significantly to mixing at the boundary. We can expect 
that these scales are distributed spatially; the smallest scales are 
found deeper in the overshooting region. This physical picture of 
the scales involved in overshooting is the basis for the interpre¬ 
tation of radiative effects outlined in the next section. 


3.3. Regimes of overshooting depending on the importance 
of radiative effects 

A turbulent flow is characterized by a wide range of eddies with 
velocity scale u and length scale /. Although the proper definition 
of a “turbulent Peclet number” is not trivial, we can consider Eq. 
(7) as a qualitative definition of a scale-dependent Peclet num¬ 
ber. A critical scale is the scale Id that corresponds to Pe ~ 1. 
Below this scale, eddies are strongly affected by thermal diffu¬ 
sion, above this scale, eddies evolve adiabatically. The assump¬ 
tion that the transition occurs exactly at Pe = 1 is a simplifi¬ 
cation, but it is sufficient for the picture presented in this work. 
Given the range of scales [/min> Wx] that contribute to overshoot¬ 
ing, we define three different regimes of overshooting depending 
on the location of the transition Pe ~ 1, as summarized in Fig. 1 
and detailed below. 


3.3.1. Case 1; /max < Id - Diffusion-dominated regime 

In this case, all the scales participating in the overshooting pro¬ 
cess have Pe < 1 and are strongly affected by thermal diffusion. 
Turbulent eddies are only able to mix composition, without af¬ 
fecting the entropy structure significantly (see Zahn 1991). This 
is the spirit of the “chemical mixing” parametrization presented 
in Sect. 2.1. This parametrization was first proposed in Freytag 
et al. (1996). Freytag (1995) presents radial profiles of the Peclet 
number for the A-type star and white dwarf models published in 
Freytag et al. (1996) (see their Figs. 11 & 12). They show that 
Pe < 10 across the convective and overshooting region, which is 
consistent with our picture. 
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3.3.2. Case 2: Zmin < Id < ^max - Penetrative regime 

In this case, there is a transition within the boundary layer: scales 
that have Pe > 1 mix entropy and composition, whereas scales 
that have Pe < 1 are only able to mix composition. This struc¬ 
tures the boundary layer in two parts: a nearly adiabatic, well- 
mixed sublayer, and a “diffusion-dominated” sublayer where 
only the composition is mixed (see Zahn 1991). This picture is 
supported by the structure of the boundary layer at the bottom 
of the convective envelope of the red giant model described in 
Viallet et al. (2013). Studying the mean-held budget equation 
for internal energy, the authors hnd that the boundary layer is 
in a thermal balance in which the heat advected by entropy-rich 
material into the layer is counterbalanced by the cooling by ra¬ 
diative diffusion. This results locally in a “super-stellar" radiative 
flux (Li-ad > T*). 

We now show that the boundary layer is indeed characterized 
by a transition in the turbulent Peclet number, consistent with the 
picture outlined above. We hrst introduce the so-called Ozmidov 
scale: 



where e is the turbulent dissipation rate of kinetic energy, N the 
Brunt-Vaisala frequency, and Lg is the length scale of the small¬ 
est eddies that are affected by buoyancy (e.g., Ozmidov 1965; 
Smyth & Mourn 2000). We then introduce the length scale re¬ 
sulting from equating the thermal diffusion time scale to the 
buoyancy time scale: 


l} 



(9) 


where can be interpreted as the scale below which ther¬ 
mal diffusion “erases” the effect of buoyancy. Inspired by the 
literature on turbulent mixing (Ozmidov 1965; Osborn 1980; 
Brethouwer & Lindborg 2009), we take as a proxy for the 
turbulent transport of heat. One therefore has 


( 




( 10 ) 


Figure 2 shows the profiles of Lg, L^, and Pe in the boundary 
layer of the red giant model of Viallet et al. (2013). There is 
indeed a transition in the Peclet number in the boundary layer 
from ~ 10^ to ~ 10 ', where the length scales Lg and cross 
each other. 

The existence of such a transition is key because it allows a 
balance to be reached: overshooting extends the size of the adi¬ 
abatic region until radiative effects inhibit the process^. We are 
currently not able to predict the depth at which this boundary 
layer develops, since it would require a model for the turbulent 
dissipation of plumes. Theoretical predictions for the extent of 
the penetration region exist (e.g., Schmitt etal. 1984;Zahn 1991; 
Hurlburt et al. 1994; Rempel 2004), but they are based on the 
rather laminar picture of plumes decelerating due to buoyancy 
braking alone. For the moment, one can rely on the parametriza- 
tion for “penetrative convection” described in Sect. 2.2. This 

^ This local thermal balance is reached quite rapidly, namely on a few 
dynamical time scales. Naturally, as the stellar structure evolves glob¬ 
ally on a Kelvin-Helmholtz (thermal) time scale, the size of the over¬ 
shooting region will evolve. 


parametrization does not include the diffusion-dominated layer 
in which composition is mixed. Zahn (1991) argues that this 
layer is very thin (see, however, Rempel 2004). In the red-giant 
model of Viallet et al. (2013), it occupies roughly one third of 
the boundary layer. This could result from the limited resolution 
inherent to numerical simulations of stellar interiors and requires 
further investigation. 

Finally, it should be mentioned that Brummell et al. (2002) 
do not find any evidence of penetrative convection in their 3D 
simulations. This could be due to their somewhat artificial setup, 
in which the profile of thermal conductivity is hxed and/or to 
the fact that their direct numerical simulations probe a different 
dynamical regime given the values of the dimensionless numbers 
(Prandtl, Reynolds, Rayleigh numbers) that they could afford in 
their investigation (Rempel 2004). 

3.3.3. Case 3: Id < Km - Entrainment regime 

In this case, all relevant scales have Pe » 1, and the overshoot¬ 
ing process can be considered to be completely adiabatic. The 
process of mass entrainment results in a growth of the well- 
mixed (entropy-Hcomposition) region. The key difference with 
the previous case is that there is nothing that can counterbalance 
the process: entrainment proceeds as long as energy is injected 
into the system. 

This regime is observed by Meakin & Arnett (2007) in their 
calculations of the convection in the oxygen-burning shell of 
a massive star. The reason is that late phases of stellar evolu¬ 
tion are driven by neutrino cooling, which acts on a time scale 
shorter than photon cooling (Arnett 1996). As a result, radiative 
diffusion becomes negligible and does not affect overshooting. 
Although their boundary layers have finite thickness with a non¬ 
trivial internal structure, the authors have shown that their data 
is described well by the parametrization of “bulk” entrainment 
presented in Sect. 2.3. 

Turbulent entrainment is also observed in hydrodynamical 
simulations of the He-flash by Mocak et al. (2009, 2010, 2011). 
This specihc phase of stellar evolution is strongly out of ther¬ 
mal equilibrium owing to the ignition of helium in a degenerate 
environment. The large release of nuclear energy drives convec¬ 
tion on a time scale that is much shorter than the radiative time 
scale. Radiative diffusion becomes negligible and turbulent en¬ 
trainment is obtained, as in the previous case. 


4. Conclusion 

Based on previous theoretical and numerical work, this paper 
presents a simple physical picture of the effects of radiation on 
overshooting. Other non-adiabatic effects exist in stars (nuclear 
burning, neutrino losses), but from current simulations we do 
not see any evidence of a direct role for these effects on over¬ 
shooting. We show that the three existing parametrizations for 
overshooting correspond to different regimes, depending on the 
importance of radiative effects. As a result, we suggest that these 
three different parametrizations should coexist in ID stellar evo¬ 
lution codes and be applied consistently depending on the lo¬ 
cal physical conditions. Ideally, the selection of the adequate 
parametrization should be based on the value of the Peclet num¬ 
ber in the boundary layer. However, one faces two difficulties: 

1. How do we estimate the Peclet number? According to our 
work, the quantity that allows distinguishing between the 
different thermal regimes is given by Eq. (10). However, 
this quantity cannot be computed in the framework of MET, 
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which has, by construction, e - 0 (no flow) when > 0 . 
The Peclet number within the convective region can be esti¬ 
mated easily with the formula 


T 

where Z)mlt = jWmlt^mlt is the usual diffusion coefficient 
computed from MLT. In a first approach, this quantity could 
be used to extrapolate the value of the Peclet number in the 
boundary layer. 

2. Given a definition of the Peclet number at the boundary, what 
are the range of values that define each regime? The simple 
picture developed in this work allows us to distinguish be¬ 
tween three thermal regimes of overshooting, depending on 
whether the effect of radiation are dominant, marginal, or 
negligible. However, it does not make any prediction about 
the transition between these regimes, whether in terms of the 
values of the Peclet number at with such transitions occur or 
in terms of how sharp such transitions are. 

Therefore, we are not yet able to give a practical criterion 
that would allow automatic selection of the appropriate descrip¬ 
tion during a stellar evolution calculation. However, as a rule of 
thumb, we suggest 

1 . applying “chemical mixing” near the surface of stars, where 
radiative effects become important (inefficient convection); 

2 . applying “penetrative convection” in the deep interior of 
stars (efficient convection) for the phases which are photon- 
cooled (core-overshooting on the main-sequence, convective 
envelope undershooting, etc.); 

3. applying “turbulent entrainment” in the deep interior when 
the evolution is driven by neutrino losses (late stages of stel¬ 
lar evolution) or in phases that are in strong thermal imbal¬ 
ance (e.g., He-flash in low-mass stars). 

It should be stressed that we do not claim that the curi'ent 
parametrizations describe each regime accurately. In particular, 
the calibration of free parameters is still required. However, a 
good starting point is to use the relevant parametrization for each 
regime. For instance, it is very likely incorrect to apply turbulent 
entrainment to describe convective cores on the main sequence 
as done in Staritsin (2013), since in this phase the time scales are 
long enough for radiative diffusion to affect overshooting. Like¬ 
wise, considering our current understanding, it is inconsistent to 
use the “chemical mixing” description deep inside the stellar in¬ 
terior; this may contribute to the errors found by Schindler et al. 
(2014). 

Failure to reproduce observational trends should then point 
to flaws in the parametrizations of these different regimes, such 
as those due to missing physics. For instance, the effects of com¬ 
position are curi'ently taken into account very poorly, mainly be¬ 
cause of the (arbitrary) choice of using the Schwarzschild or the 
Ledoux criteria and the ad hoc use of thermohaline mixing. In 
fact, like radiation, composition probably affects the regime of 
overshooting and could drastically change the simple picture of 
the three regimes defined in this work. Likewise, rotation and 
magnetic held could also affect the nature and the efficiency of 
mixing at convective boundaries. 

Here as well, insight from multidimensional simulations 
should help us to understand which are the relevant parameters 
that characterize these processes and how the parameter space 
is split into different regimes of overshooting. This work was a 


first step toward this goal, and it was focused on thermal diffu¬ 
sion. Applying this strategy to other physical processes should 
bring us closer to a physically consistent description of convec¬ 
tive boundary mixing in ID stellar evolution codes. 
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